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Abstract 

We propose an algorithm for computing the proximity operator of a sum of composite convex func- 
tions in Hilbert spaces and investigate its asymptotic behavior. Applications to best approximation 
and image recovery are described. 

Keywords: Best approximation, convex optimization, duality, image recovery, proximity 
operator, proximal splitting algorithm 



1. Introduction 

in ! 

^r) ', Let % be a real Hilbert space with scalar product (• | •) and associated norm || • ||. The best 

approximation to a point z £ H from a nonempty closed convex set C C % is the point Pqz £ C 
that satisfies ||Pc£ — ^|| = rnhixec \\x — z\\. The induced best approximation operator Pq'- H —¥ C, 
also called the projector onto C, plays a central role in several branches of applied mathematics 



lOp . If we designate by lc the indicator function of C, i.e., 

if x G C; 

if. .r i—' < 




(1.1) 



then Pqz is the solution to the minimization problem 

minimize lc{x) H — ||x — z\\ . (1-2) 

x£+t 2 
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Now let Tq(T-L) be the class of lower semicontinuous convex functions f:H—>] — oo, +00] such that 



dom/ = {x G % | f(x) < +00} ^ 0. In 13] Moreau observed that, for every function / G ro(%), 
the proximal minimization problem 

minimize f(x)-{ — \\x — z\\ (1-3) 

xeH 2 

possesses a unique solution, which he denoted by proxj z. The resulting proximity operator 
prox 1 : 7i — > 71 therefore extends the notion of a best approximation operator for a convex set. 
This fruitful concept has become a central tool in mechanics, variational analysis, optimization, 
and signal processing, e.g., [ll.l7l.ll6l]. 

Though in certain simple cases closed-form expressions are available [7j, |g, Il4 ], computing 



proxj- z in numerical applications is a challenging task. The objective of this paper is to propose 
a splitting algorithm to compute proximity operators in the case when / can be decomposed as a 
sum of composite functions. 

Problem 1.1 Let z G H and let (wj)i<j< m be reals in ]0, 1] such that X^i w * = 1- F° r ever y 
i G {1, . . . , m}, let (Qi, \\ ■ ||gj be a real Hilbert space, let ri £ Qi, let gi G ro(^j), and let Lj : % — >■ Qi 
be a nonzero bounded linear operator. The problem is to 

m 

minimize > Uigi(LiX — r^) -\ — ||x — z\\ . (1-4) 

i=l 

The underlying practical assumption we make is that the proximity operators (prox q .)i<i< m 
are implementable (to within some quantifiable error). We are therefore aiming at devising an 
algorithm that uses these operators separately. Let us note that such splitting algorithms are 
already available to solve Problem 11.11 under certain restrictions. 

A) Suppose that Q\ = H, that L\ = Id , that the functions (gi)2<i<m are differentiable ev- 
erywhere with a Lipschitz continuous gradient, and that r^ = 0. Then ([1.4j> reduces 
to the minimization of the sum of f\ = g\ G Toffl) and of the smooth function fo 
s ^A=2 bJ i9 i ° Li + || ■ — -2|| 2 /2, and it can be solved by the forward-backward algorithm [81. 1 IE 



B) The methods proposed in [J] address the case when, for every i G {l,...,m}, Qi = 7i, 
Li = Id , and rj = 0. 

C) The method proposed in [5] addresses the case when m = 2, Q\ = %, and L\ = Id , and 
n = 0. 

The restrictions imposed in A) are quite stringent since many problems involve at least two non- 
differentiable potentials. Let us also observe that since, in general, there is no explicit expression 
for prox 9 . oLj in terms of prox fl . and Lj, Problem 11.11 cannot be reduced to the setting described 



in B). On the other hand, using a product space reformulation, we shall show that the setting 
described in C) can be exploited to solve Problem 11,11 using only approximate implementations of 
the operators (prox 9 .)i<j< m . Our algorithm is introduced in Section [21 where we also establish 
its convergence properties. In Section [3l our results are applied to best approximation and image 
recovery problems. 

Our notation is standard. ¥> (H,G) is the space of bounded linear operators from H to a real 
Hilbert space Q. The adjoint of L € 23 (%, Q) is denoted by L*. The conjugate of / E Tq(T-L) is the 
function /* G Tq(T-L) defined by /*: u i-> sup a , eW ((2; | u) — /(x)). The projector onto a nonempty 
closed convex set C C % is denoted by Pq- The strong relative interior of a convex set C C % is 



sriC = {x € C | cone(C — x) = span (C — x)}, 



where coneC = M {Ax | x € C}, (1.5) 



A>0 



and the relative interior of C is riC = {x € C cone(C — x) = span (C — x)}. We have int C C 
sri C C ri C C C and, if % is finite-dimensional, ri C = sri C. For background on convex analysis, 
see hi 



2. Main result 

To solve Problem 11.11 we propose the following algorithm. Its main features are that each 
function gi is activated individually by means of its proximity operator, and that the proximity 
operators can be evaluated simultaneously. It is important to stress that the functions {gi)\<i< m 
and the operators {Li)i<i< m are used at separate steps in the algorithm, which is thus fully 
decomposed. In addition, an error aj iTl is tolerated in the evaluation of the ith proximity operator 
at iteration n. 

Algorithm 2.1 For every i G {1, . . . ,m}, let (ai ;n ) nS N be a sequence in Q\, 
Initialization 



P 



maxi<j< m \\Li 



-2 



e € ]0,min{l, p}[ 

For i = 1 , . . . , m 
[ v i>0 E Qi 

For n = 0, 1, . . . 

Em j * 

7„ G [e, 2/) - e] 

An € [e, 1] 

For i - 1 , . . . , m 

!',.,,+ 1 = Vj, n + A n ( prox 7n;/ 



(2.1) 



( prox 7nS * (v i)n + 7„(LjX n - r,)) + a ijn - Uj jn J , 



Note that an alternative implementation of ()2.ip can be obtained via Moreau's decomposition 
formula in a real Hilbert space Q [8J, Lemma 2.10] 

(V5er (£))(V 7 €]0,+oo[)(Vi;e£) prox 7g *t; = U -7prox 7 -i 9 (7- 1 T;). (2.2) 

We now describe the asymptotic behavior of Algorithm 12.11 
Theorem 2.2 Suppose that 

(n)i<i< m G sri{(Lj2;-yj)i<j< m | x G U,(yi)i<i< m G X^Lidom^} 



(2.3) 



and that 



(Vi e {l,...,m}) ^ IKnlb < 

nGN 



-CO. 



(2.4) 



Furthermore, let (x n ) n ^n, (v i, n )neN> • • • > (%,n)neN ^ e sequences generated by Alaorithm \2.1\ Then 
Problem \l.l\ possesses a unique solution x and the following hold. 

(i) For every i G {1, . . . ,m}, (vi, n ) n eN converges weakly to a point V{ G £/«. Moreover, (uj)i<i< m 
is a solution to the minimization problem 



1 

minimize — 

Bl6ffl »m6Sm 2 



y^ ^iL*^ 



t=i 



+^ wi (^*(^) + fa i ri ))' 



(2.5) 



i=l 



and x = z - Xl^Li WjI/*Uj. 

(ii) (x n ) n£ N converges strongly to x. 

Proof. Set f:W—} ]— oo,+co] : x t-> Y^i^idii^iX ~~ r *)- The assumptions imply that, for every 
i G {1, . . . , m}, the function x i— » gi{L%x — r,) is convex and lower semicontinuous. Hence, / is 
likewise. On the other hand, it follows from (|2.3p that 



(r"i)l<i<m G {(ijX-yj)i<j< m | x G V.,(yi)i<i< m G X^Lidom^} 



(2.6) 



and, therefore, that dom/ ^ 0. Thus, / G To(^) and, as seen in (jl.3p . Problem 11.11 possesses a 
unique solution, namely x = prox<- z. 

Now let T-L be the real Hilbert space obtained by endowing the Cartesian product "H m with the 
scalar product (• | -) H : (x,y) i-> XXi^j^i | yi), where a; = (xj)i<j< m and y = (yi)i<i< m denote 
generic elements in T-L. The associated norm is 



\n- x^ 



Ell 112 
UJi\\Xi\\ 2 . 

\ i=l 



(2.7) 



Likewise, let Q denote the real Hilbert space obtained by endowing the Cartesian product Q\ x 
• • • x Q m with the scalar product and the associated norm respectively defined by 



(• I -)g- (y>«) ^^2^i{Vi \ Zi)g. and \\-\\g-V^ 



i=\ 



\ i=l 



Define 



f = ld, where D = {(x, . . . ,x) G ~H | i£ %} 
g: Q ->-]-oo,+oo] : y h-> Y!iLi^i9i{Vi) 
L-.'H^Q-.x^f (LiXi)x<i<m 

r = (ri,...,r m ) 
z= (z,...,z). 



(2.J 



(2.9) 



Then / € r'o('H), 9 € To(^), and L € !B {H,Q). Moreover, D is a closed vector subspace of T-L 
with projector 



proxj = P D : x i-fr I y^ WjXj, . . . , ^ w;Xj 



i=l 



and 



L* : Q ^ U: v ^ (L*^)_. ,., . 
V * vl<i<m 

Note that ([23]) and fl27T]) yield 



(2.10) 



(2.11) 



(VxeH) \\Lx\\g =y^ j UJi\\L i Xi\\g. 



1=1 

m 

4 = 1 



2 II i|2 



< I max IlLjII ) > 

\l<Km / ^^ 

4 = 1 

max ||Lj|| )||aj||-w- 

KKra / 



U; Xi 



Therefore, 



_L| < max \\L{ 

Ki<m 



(2.12) 



(2.13) 



We also deduce from (|2.3[) that 

r G sri (L(dom/) — doing). 
Furthermore, in view of (|2.7p and (|2.9p . in the space T-L, (jl.4p is equivalent to 

minimize f(x) + g(Lx — r) -\ — \\x — z\fo. 
xe~H 2 

Next, we derive from [a, Proposition 3.3] that the dual problem of (|2.15p is to 
minimize f*(z — L*v) + g*{v) + (v | r) g , 

where f :«4 inf^g-^ (/*(«;) + (1/2) ||w — w||;L) is the Moreau envelope of /*. Since 
we have /* = i D ± . Hence, (|2.7p and (|2.10p yield 

m '■ 



(Vt* g tt) /* (t») = i||« - Pd-lwII^ = ^ll^ll^ = ~ 



i=l 



On the other hand, ([2TH|) and (J^ISJ) yield 

m 

(Vw G </) g*(f) = y^^iff*fa) and prox g * v = (prox s * 



Vi )i<i< m - 



i=l 



Altogether, it follows from (gUp) , ([2~TTD . (l2~T8|) . and (USD, that 

(J27L6]) is equivalent to (|23|) . 
Now define 



(Vn G N) 



Xn — \%m ■ ■ ■ t %n) 
V n = (Vl.Ti)- • • |fm,n) 



(2.14) 

(2.15) 
(2.16) 

/ = LD, 

(2.17) 

(2.18) 
(2.19) 

(2.20) 



a, 



(ai n , . . . , a 



n , . . . , ui m ,nj 



Then, in view of (JZHJ), (l2J0]l . (I2TTTD . (I2TT3D . and (I2TT8D . (l2~Tjl is a special case of the following 
routine. 

Initialization 

p=\\l\\- 2 

e G ]0,min{l, p}[ 
v £G 
For n = 0,1,... (2.21) 

s n = prox f (z - L*v n ) 
7„ G [e, 2p - e] 

An G [e, 1] 

v n +i = v n + A„ ( prox 7ng , (v n + 7 n (Lx n -r))+a n -v n ). 



Moreover, (|2.4|) implies that X^ngN ll a nlle < +00. Hence, it follows from (|2.14p and [5|, Theorem 3.7] 
that the following hold, where x is the solution to (|2.15p . 

(a) (t> n ) raG N converges weakly to a solution v to (|2.16|) and x = proxj(2 — L*v). 

(b) (x n ) n( z^ converges strongly to x. 

In view of ([22]), (jMJ, ([2Jj]), (I2TTUD . (I2TT11 . (f2TT9l) . and (12301) . items (a) and (b) provide respectively 
items 



(i) and (ii) □ 



Remark 2.3 Let us consider Problem 11.11 in the special case when (Vi G {l,...,m}) £/j = H, 
Lj = Id , and r, = 0. Then (|1.4p reduces to 



minimize 



rra 1 

^Wi£/i(x) + -||x - z\\ 2 . 



(2.22) 



8=1 



Now let us implement Algorithm 12.11 with j n = 1, A n = 1, aj jn = 0, and Vj : o = 0. The iteration 
process resulting from (|2.ip can be written as 



(2.23) 



Initialization 
xo = z 
For i = 1 , . . . , m 

For n = 0, 1, . . . 
For i = 1, . . . ,m 

Vi, n +i = prox g «(x n + v i>n ). 

X n +l = Z - YT=l u i v i,n+l- 

For every % G {1, . . . , m} and n G N, set £j ; „ = x n + i^ )n . Then (J2.23P yields 

Initialization 

x = z 

For i = 1 , . . . , m 

[ Zifl = z 
For n = 0, 1, . . . 

^n+1 = Z- YT=l U i P rox 9 * ^,n 

For i = 1, . . . , m 

Zi,n+1 = Xn+l + PrOXg* 2^. 

Next we observe that (Vn G N) X^li u i z i,n = £• Indeed, the identity is clearly satisfied for 



(2.24) 



;; 



and, for every n G N, (12.24J) yields ^^i w i^,n+i = %n+i + YaL\ ^i prox 9 » z i>n 



Y^ILi u i P rox a* z in) + YllLi u i P rox o* z i n = %• Thus, invoking (|2.2p with 7 = 1, we can rewrite 
(12321) as 

Initialization 
x = z 

For i = 1, . . . ,m 
I ^,o = z 

T7 n 1 ( 2 - 25 ) 

For n = 0, 1, . . . v ' 

x n+l = YT=1 ^i P r ° x g, Zi,n 

For % = 1, . . . ,m 

|_ -2j,n+l — 2-n+l T 2j,n P^X^ ^i,n- 

This is precisely the Dykstra-like algorithm proposed in [4j, Theorem 4.2] for computing 
proxy-™ w z (which itself extends the classical parallel Dykstra algorithm for projecting z onto 
an intersection of closed convex sets [2J, [llj). Hence, Algorithm 12.11 can be viewed as an extension 
of this algorithm, which was derived and analyzed with different techniques in [4| . 

3. Applications 

As noted in the Introduction, special cases of Problem 11.11 have already been considered in the 
literature under certain restrictions on the number m of composite functions, the complexity of 
the linear operators (Lih<i<mi and/ or the smoothness of the potentials (gi)x<i<m (one will find 
specific applications in [3, la, 13, la, |9|, LL5[ and the references therein). The proposed framework 
makes it possible to remove these restrictions simultaneously. In this section, we provide two 
illustrations. 

3.1. Best approximation from an intersection of composite convex sets 

In this section, we consider the problem of finding the best approximation Ppz to a point 
zeM from a closed convex subset D of T-i defined as an intersection of affine inverse images of 
closed convex sets. 

Problem 3.1 Let z € H and, for every i € {1, . . . , rn\, let (Qi, \\ ■ ||g.) be a real Hilbert space, let 
n G Qi, let Ci be a nonempty closed convex subset of Qi, and let ^ Li £ 15 (%,Qi). The problem 
is to 

in 

minimize \\x — z\\, where D = (| [x € H I LiX € T{ + CA. (3-1) 

X&D II 1 -! 

1=1 

In view of (jl.ip . Problem l3.1l is a special case of Problem ll.lt where (Vi £ {1, . . . , m}) gi = tc t 
and Ui = 1/m. It follows that, for every iG{1,,.. ,m} and every 7 E ]0, +00 [, prox 79 . reduces to 
the projector Pq 4 onto Cj. Hence, using (|2.2p . we can rewrite Algorithm 12. II in the following form, 
where we have set Ci >n = — 7„ Oi )7l for simplicity. 



Algorithm 3.2 For every i G {1, . . . , m}, let (cj jn ) ne N be a sequence in Gi 



Initialization 

p= (maxi<j< m ||Li||) 
e G ]0,min{l, p}[ 
For i = 1 , . . . , m 
[ v ifi G & 

For n = 0, 1, . . . 

Em j * 

7„ G [e, 2p - e] 

An G [e, 1] 

For i = 1 , . . . , m 

Ui,„ + 1 = V itn + 7„A n f LjX n -n- Pc, {ln lv i,n + LiX n - n) - Ci,n) ■ 



(3.2) 



In the light of the above, we obtain the following application of Theorem 12. 2;(ii) 
Corollary 3.3 Suppose that 



{ri)i<i<m G sri{(L ; 



X ~ Vi)l<i<r 



s H, (yi)i<i< m G X^LiCj} 



(3.3) 



nGN 11^," 



|^ < +oo. TTien every sequence (x n ) ng N generated by 



and that (Vi G {1, . . . ,m}) ^ 

Alaorithm \3.2\ converges strongly to the solution Pr>z to Problem \3.1 

3.2. Nonsmooth image recovery 

A wide range of signal and image recovery problems can be modeled as instances of Problem ll.il 
In this section, we focus on the problem of recovering an image x G 7i from p noisy measurements 



TiX + Si, 1 < i <p. 



(3.4) 



In this model, the ith measurement rj lies in a Hilbert space Gi, Pi G "B (T~L, Gi) is the data formation 
operator, and Sj G Gi is the realization of a noise process. A typical data fitting potential in such 
models is the function 



t-> \ L0igi(TiX — rj), where < $j G Fq(^j) and <?j vanishes only at 0. 



(3.5) 



i=i 



The proposed framework can handle p > 1 nondifferentiable functions (gi)i<i< p as well as the 
incorporation of additional potential functions to model prior knowledge on the original image x. 
In the illustration we provide below, the following is assumed. 



The image space is H = Hq(Q), where Q is a nonempty bounded open domain in M 2 . 

x admits a sparse decomposition in an orthonormal basis (ek)keN of H. As discussed in 
[a, [20J this property can be promoted by the "elastic net" potential x i-> ^CfceN^fc(( x I £k))i 
where (V/c € N) 4>k : ■? *- > a\£\ + @\£\ 2 , with a > and /3 > 0. More general choices of suitable 
functions (</>fc)fceN are available [6]. 

x is piecewise smooth. This property is promoted by the total variation potential tv(x) = 



f n \Vx(uj)\2duj, where | • I2 denotes the Euclidean norm on M. 2 17. | 



Upon setting gi = || ■ ||ft in (|3.5p . these considerations lead us to the following formulation (see 
[a, Example 2.10] for more general nonsmooth potentials). 

Problem 3.4 Let % = Hg(f2), where O C K 2 is nonempty, bounded, and open, let {wi)i<i< p +2 
be reals in ]0, 1] such that X^f=i w « = 1j anc ^ ^ ( e fc)fcgN ^ e an orthonormal basis ofH. For every 
i € {1, . . . ,p}, let 7^ Tj G 23 (%, C/j), where (Qi, \\ ■ ||gj is a rea/ Hilbert space, and let r% G <5«. T7te 
problem is to 

y^,Wi\\TjX - n^ +y~^ (u p+1 \{x I ejfc)| + ~|(cc I e fe )| 2 J +w p+2 tv(j;). (3.6) 



minimize 



It follows from Parseval's identity that Problem 13.41 is a special case of Problem 11.11 in % = 
Hq(O) with m = p + 2, z = 0, and 

9i=\\- ||ft and Lj = Tj, if 1 < i < p; 

Gp+i = £ 2 (N), g p+ \ = || • ||^i, r p+1 = 0, and L p+1 : x i-+ ((x | e fc )) fceN ; (3-7) 

,(?p+2 = L 2 (fi) L 2 (fi), 5 P+2 : U ^ In l2/( w )Mw, r p+2 = 0, and L p+2 = V. 

To implement Algorithm l2.il it suffices to note that L* +1 : (v^keN ^ YlkeN^^k and L* +2 = 
— div, and to specify the proximity operators of the functions (jg*)i<i<. m , where 7 G ]0, +oo[. First, 
let i € {1, . . . ,p}. Then g^ = \\ ■ \\g. and therefore g* = Ib v where L>i is the closed unit ball of Qi. 
Hence prox 7g . = Pb. . Next, it follows from (|2.2p and [8j, Example 2.20] that prox 7ff * : (^OfceN ^ 



(-P[-i,i]6c)fcgN- Finally, since g p +2 is the support function of the set pi 

K = {yeg p+2 I |y|2 < 1 a.e.}, (3.8) 

5* +2 = Ik and therefore prox 7g * = Pk, which is straightforward to compute. Altogether, as 
||-L p+ i|| = 1 and ||L P +2|| < 1> Algorithm 12.11 assumes the following form (since all the proximity 
operators can be implemented with simple projections, we dispense with the errors terms). 



10 



Algorithm 3.5 

Initialization 

p= (max{l,||ri||,...,||T p ||})~ 
e G ]0,min{l,p}[ 
For i = 1, . . . ,p 
[ fj,0 G £i 

Vp+1,0 = (vk,o)keN G ^ 2 (N) 

V2,o£L 2 (l])0L 2 (fi) 

For n = 0, 1, . . . 

^n = Z - Yfi=l UiT*Vi^ n - UJ p+1 J2keN v k,nZk + ^p+2 div V p+2 ,n 

7ne[e,2p-e] (3-9) 

An G [e, 1] 

For 2 = 1, . . . ,p 

, / ^,n+7n(^n-^i) \ 

fj,n+l — ^j,n + *n I f1 ,, , \\\ ^ v i,n) 

„ , ,_ m . / ^,n + 7n(^n I e fc ) 

For every k G N, i>fc, n +l = i/ fc n + A n t— ■ ■ : -p- - v kn 

\max{l, \Vk,n + Jn{Xn I e k )\} 

For almost every u G f2, 

"p+2,n+lW = fp+ 2 ,n (w) + A n r— -^— - -^— - - V p+2 ,n{0J) ■ 

Vmaxjl, |u p+2 ,nM +7nVx n (a;)|2} F / 

Let us establish the main convergence property of this algorithm. 

Corollary 3.6 Every sequence (a%)neN generated by Algorithm \3.5\ converges strongly to the so- 
lution to Problem\3.4\ 



Proof. In view of the above discussion and of Theorem I2.2;(ii)[ it remains to check that (12.3|) is 
satisfied. Set S = {(LiX - yi)i<i< m \ % G H, (yi)i<i< m G X™idom#i}. We have domgj = Qi for 
every i G {1, . . . ,p}, domg p+ i = ^ 1 (N), and domg p+ 2 = L 2 (fi) © L 2 (fi). Consequently, 

S = UTix -yi,.. .,T p x - y p , ({x \ e k ) - rj k )ken, Vx - y p+2 

xeU, {yi)x<i< p G X p i=1 Gi, (%) fee N G ^(N), y P+2 G L 2 (ft) © L 2 (ft)} 
= (Xf =1 &) x ^( N ) x (L 2 (fi) ©L 2 (fi)) 

= X£i&. (3.10) 

Hence, we trivially have (n, . . . , r p , 0, 0) G sri S. D 

11 



Let us emphasize that a novelty of the above variational framework is to perform total variation 
image recovery in the presence of several nondifferentiable composite terms, with guaranteed strong 
convergence to the solution to the problem, and with elementary steps in the form of simple 
projections. The finite-dimensional version of the algorithm can easily be obtained by discretizing 
the operators V and div as in [3J (see also [5j, Section 4.4] for variants of the total variation 
potential). 
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